// ModifiedGriewankFunction.cpp: implementation of the ModifiedGriewankFunction class.
//
//      MGr = ((x1^2 + x2^2)/4000) - Cos(x1)Cos(x2/Sqrt(Period)) + 1
//
// Used to hopefully help us identify which masses to use!
//////////////////////////////////////////////////////////////////////////////////////
#include <iomanip>
#include <fstream>
#include <string>
#include "modifiedgriewank_function.h"
#include "DataHolder.h"



//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

ModGriewankFunction::ModGriewankFunction(int i_NInput, double Period)
{
  NInput = i_NInput;
  m_Period = Period;
  nCall = 0;

  newA(m_Derivative,NInput);

  Result.Size(NInput);
}


ModGriewankFunction::~ModGriewankFunction()
{
  deleteA(m_Derivative);
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void ModGriewankFunction::checkBoundary(vector<double> &Input)
{
  for (unsigned int i = 0; i < Input.size(); i++)
  {
    // periodic boundary
    if ( Input[i] >= 100.0 )
    {
      Input[i] = -100.0 + fmod(Input[i],100.0);
    }
    else if ( Input[i] <= -100.0 )
    {
      Input[i] = 100.0 + fmod(Input[i],100.0);
    }
  }
}

// takes a point and if it is beyong boundary move the point
// into the boundary box
void ModGriewankFunction::checkBoundary(double *Input)
{
  for (int i = 1; i <= NInput; i++)
  {
    // periodic boundary
    if ( Input[i] >= 100.0 )
    {
      Input[i] = -100.0 + fmod(Input[i],100.0);
    }
    else if ( Input[i] <= -100.0 )
    {
      Input[i] = 100.0 + fmod(Input[i],100.0);
    }
  }
}

//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double ModGriewankFunction::Value(double *Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  if(NInput != 2)
  {
    cout << "Sorry - you have tried to use this function \nwith the wrong number of variables"  
      << "\nthe only number is two" << endl;
    exit(1);
  }

  double danswer = 0.0;
  double answer = 0.0;

  for(int i = 1; i<=NInput; i++)
  {
    answer += (Input[i]*Input[i])/4000;
  }

  danswer += cos(Input[1]);
  danswer *= cos(Input[2]/m_Period);

  answer -= danswer;
  answer += 1;

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// calculates/builds the function value
double ModGriewankFunction::Value(vector<double> &Input)
{
  // first check that point is not beyong any boundaries
  checkBoundary(Input);

  if(NInput != 2)
  {
    cout << "Sorry - you have tried to use this function \nwith the wrong number of variables"  
      << "\nthe only number is two" << endl;
    exit(1);
  }

  double danswer = 0.0;
  double answer = 0.0;

  for(int i = 1; i<=NInput; i++)
  {
    answer += (Input[i]*Input[i])/4000;
  }

  danswer += cos(Input[1]);
  danswer *= cos(Input[2]/sqrt(double (m_Period)));

  answer -= danswer;
  answer += 1;

  Result.Input(answer, Input);

  nCall++;

  return (answer); 
}

// calculates/builds the derivative (used for CGM and HMC)
void ModGriewankFunction::BuildDerivative(double *Input, double *Derivative)
{
  for(int j = 1; j<=NInput; j++)
  {
    Derivative[j] = 0.0;
  }

  Derivative[1] += (Input[1]/2000); 
  Derivative[1] += sin(Input[1])*cos(Input[2]/m_Period);

  Derivative[2] += Input[2]/2000;
  Derivative[2] += ((1.0/m_Period)*sin(Input[2]/m_Period))*cos(Input[1]);

  for(int k = 1; k<=NInput; k++)
  {
    m_Derivative[k] = Derivative[k];
  }
}

// used to test if the function is working properly
void ModGriewankFunction::test_function()
{
  ofstream out("output/test_of_modifiedgriewank.txt");
  
  out << "NUMBER_OF_DIMENSIONS: " << NInput << endl 
    << "      Subscript[x,k]      |     Value    |         Derivative[k]" << endl << endl;

  int NRun = 3;
  double *mInput;
  double *mDeriv;
  newA(mInput, NInput);
  newA(mDeriv, NInput);

  for(int run=1; run<=NRun; run++)
  {
    for(int n=1; n<=NInput; n++)
    {
      mInput[n]=(200.0*ran1())-100.0;
    }

    for(int k=1; k<=NInput; k++)
    {
      out << setw(11) << mInput[k];
    }
      // d[1] == 0, doesn't give derivative = 0
      // as d[1] is never 0.
      out << setw(16) << Value(mInput);

      BuildDerivative(mInput, mDeriv);

    for(int l=1; l<=NInput; l++)
    {
      out << setw(16) << m_Derivative[l];
    }

    out << "\n\n";
    
  }

  deleteA(mInput);
  deleteA(mDeriv);
  out.close();
}

// used to print the solution (this is the old version, the new version is below)
void ModGriewankFunction::print_sol(double *solution, double funcVal, char filename[])
{
  ofstream out(filename);	

  out << "Number of dimensions: " << NInput << endl;
  out << "Final ModGriewank value = " << funcVal << endl;
  out << "Solution coordinates:" << endl;

  for (int i = 1; i <= NInput; i++)
  {
    out << solution[i] << "  ";
  }
 
  out << endl;

  out.close();
}

// used to print the solution to a file
void ModGriewankFunction::print_sol(char filename[])
{
  ofstream out(filename);	

  // getBestValue() returns the best value stored in the tracker
  out << "Number of dimensions: " << NInput << endl;
  out << "Final ModGriewank value = " << Result.getBestValue() << endl;
  out << "Solution coordinates:" << endl;

  // prints out the solution coordinates to a file
  Result.printSolution(out);

  out.close();
}


// distance from know global minimum which is located at the origin
double ModGriewankFunction::distanceFromGlobalMinimum(double *solution)
{
  double toReturn = 0.0;

  for (int i = 1; i <= NInput; i++)
  {
    toReturn += solution[i]*solution[i];
  }

  double toReturnS = sqrt(toReturn);
  return toReturnS;
}